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ABSTRACT 

We report results of a large set of iV-body calculations aimed to study the evolution 
of multi-mass star clusters in external tidal fields. Our clusters start with the same 
initial mass-functions, but varying particle numbers, orbital types and density profiles. 
Our main focus is to study how the stellar mass-function and other cluster parameters 
change under the combined influence of stellar evolution, two-body relaxation and the 
external tidal field. 

We find that the lifetimes of star clusters moving on similar orbits scale as T ~ 
T^H^ where Thh is the relaxation time, and the exponent x depends on the initial 
concentration of the cluster and is around x ss 0.75. The scaling law does not change 
significantly if one goes from circular orbits to eccentric ones. From the results for 
the lifetimes, we predict that between 53% to 67% of all galactic globular clusters 
will be destroyed within the next Hubble time. Low-mass stars are preferentially lost 
and the depletion is strong enough to turn initially increasing mass-functions into 
mass-functions which decrease towards the low-mass end. The details of this depletion 
are insensitive to the starting condition of the cluster and can be characterised as 
a function of a single variable, as e.g. the fraction of time spent until total cluster 
dissolution. 

The preferential depletion of low-mass stars from star clusters leads to a decrease 
of their mass-to-light ratios except for a short period close to final dissolution, when the 
mass-fraction in form of compact remnants starts to dominate. The fraction of compact 
remnants is increasing throughout the evolution. They are more strongly concentrated 
towards the cluster cores than main-sequence stars and their mass-fraction in the 
center can reach 95% or more around and after core-collapse. 

For a sample of galactic globular clusters with well observed parameters, we find a 
correlation between the observed slope of the mass-function and the lifetimes predicted 
by us. It seems possible that galactic globular clusters started with a mass-function 
similar to what one observes for the average mass-function of the galactic disc and 
bulge. 

Key words: methods: methods: numerical, stellar dynamics - globular clusters: gen¬ 
eral. 


1 INTRODUCTION 

The initial mass function (IMF) of stars is an important 
quantity for astronomy since it influences a wide variety of 
astrophysical processes, like e.g. the chemical evolution of 
galaxies, the dynamical evolution of star clusters and the 
determination of the absolute star formation rate. It is also 
important for theory, since any theory of star formation has 
to be able to reproduce it. 

One of the best environments to study the mass- 
function of stars are globular clusters, since each globular 
cluster contains a relatively large sample of stars with sim¬ 
ilar distances, chemical composition and ages. In addition, 
the fraction of binary stars is relatively low in globular clus¬ 


ters, of order 15% to 20% (Albrow et al. 2001, Elson et 
al. 1998), and is thought to effect the determination of the 
mass-function only in the cluster centers (Rubenstein & Bai- 
lyn 1999). 

Thanks to the availability of the Hubble Space Tele¬ 
scope and large ground-based telescopes, it is nowadays pos¬ 
sible to study the luminosity-function of stars in globular 
clusters from the tip of the main-sequence down to the hy¬ 
drogen burning limit and observations are have been car¬ 
ried out for about a dozen galactic globular clusters (see 
for example Paresce & de March! 2000, Piotto & Zoccali 
1999 and references therein). Although some uncertainty in 
transforming the observed luminosity function into a mass- 
function remains, it seems that most mass-functions can be 
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fitted by power-laws with slopes comparable to or some¬ 
what flatter then what one finds for the mass-function of 
low-mass stars in the galactic disc and bulge. There are indi¬ 
cations that not all clusters exhibit the same mass-function, 
although the reasons for the variation seen among different 
clusters are much less clear. Differences might arise due to 
differences in the metallicities (De Marchi & Paresce 1995) 
or as a result of the dynamical evolution of the clusters (Pi- 
otto & Zoccali 1999). 

However, in order to extract the scientifically much 
more interesting initial mass-function of stars from the ob¬ 
served mass-function, one has to know how the IMF has 
changed due to the dynamical evolution of the clusters. Fur¬ 
thermore, observations of stars in globular clusters usually 
cover only a limited range in radius, so the effect of mass- 
segregation of stars has to be understood in order to deduce 
the global mass-function from the observations. Although it 
is possible to correct for the effect of mass-segregation by 
assuming equipartition and fitting multi-mass King-Michie 
models to the observed luminosity profile of a cluster (Gunn 
& Griffin 1979), these question are probably best answered 
by dynamical simulations. 

The pioneering study in this field was made by Cher- 
noff & Weinberg (1990), who followed the evolution of multi¬ 
mass clusters surrounded by a tidal cut-off and evolving un¬ 
der the combined influence of two-body relaxation and mass- 
loss from stellar evolution by means of an Fokker-Planck 
code. They found that weakly concentrated clusters with 
a flat mass-function are easily destroyed due to cluster ex¬ 
pansion and tidal mass-loss resulting from the mass-loss by 
stellar evolution. They also presented a number of obser¬ 
vational parameters for clusters which survived up to the 
present time. Their results were later extended by Fukushige 
& Heggie (1995), Portegies Zwart et al. (1998), Giersz (2001) 
and Joshi et al. (2001). These studies generally confirmed 
the results of Chernoff & Weinberg for the cluster life¬ 
times, but showed that an isotropic Fokker-Planck code can 
give lifetimes significantly shorter than the results obtained 
by other methods for certain cases. Takahashi & Portegies 
Zwart (1998, 2000) solved this discrepancy by showing that 
an anisotropic Fokker-Planck code with a better treatment 
of the external tidal field gives good agreement with the N- 
body method. In their latter paper they conducted an ex¬ 
tensive survey of initial cluster conditions and showed that 
clusters can lose a large fraction of their mass in a short time 
due to stellar evolution and still remain bound. Spurzem & 
Aarseth (1996) also found good agreement between various 
computational methods for the pre and post-collapse evolu¬ 
tion of a cluster containing 10000 equal-mass stars. 

Vesperini & Heggie (1997) have performed a large set of 
A-body simulations of the evolution of star clusters evolving 
under the combined influence of two-body relaxation, stel¬ 
lar evolution and disc-shocking. Their clusters started with 
power-law mass-functions and they could establish a rela¬ 
tion between the changes of the mass-function and the mass 
which is already lost from the clusters. They also showed 
that the slope of the mass-function near the half-mass ra¬ 
dius will always stay close to the slope of the global mass- 
function. This result was later also obtained by Takahashi & 
Lee (2000) from Fokker-Planck simulations. Small-N simula¬ 
tions of multi-mass clusters in tidal fields were also carried 
out by de la Fuente Marcos (1995, 1996). His simulations 


showed that multi-mass clusters evolve considerably faster 
than single-mass clusters and that the form of the initial 
mass-function and the fraction of binary stars influences the 
lifetime of a star cluster. 

Globular clusters evolve due to a number of internal and 
external processes, as for example mass-loss due to stellar 
evolution, core-collapse and cluster re-expansion as a result 
of two-body relaxation and external tidal shocks from the 
parent galaxy. The realistic incorporation of all these effects 
into a dynamical simulation has proven to be a substantial 
challenge for both hardware and software since the different 
dissolution mechanisms act simultaneously and the range of 
timescales involved in the numerical integration of the prob¬ 
lem is rather large. So far, realistic simulations of globular 
clusters were hampered by the fact that A-body simulations 
with particle numbers realistic for globular clusters were not 
feasible, while at the same time it was difhcult to include all 
relevant physical effects into a Fokker-Planck or Monte Carlo 
code with the necessary realism. Furthermore, as was shown 
by the ’Collaborative Experiment’ (Heggie et al. 1998) and 
Baumgardt (2001), the scaling of A-body simulations from 
small-A clusters up to the globular cluster regime is not 
without risk, especially if many processes acting on differ¬ 
ent timescales and depending on the number of cluster stars 
in differing ways, are involved. 

With the completion of the GRAPE6 special purpose 
computers at Tokyo University (Makino 2002), the situation 
has however changed dramatically. The GRAPE6 boards 
allow the time-consuming calculation of the mutual grav¬ 
itational forces to be done in hardware, thereby making it 
possible to simulate the evolution of star clusters containing 
up to 3 — 5 • 10® stars with A-body programs like NBODY4 
(Aarseth 1999) or KIRA (Portegies Zwart et al. 1999) within 
a month of computing time. The number of stars feasible to 
simulate is thus comparable to what is found in small glob¬ 
ular clusters and the scaling of the results to real globular 
clusters, if necessary at all, involves only very small factors, 
and is therefore less risky. 

In this paper we present results of a large set of A- 
body calculations for the dynamical evolution of globular 
clusters done on GRAPE6. Our simulations include stellar 
mass-loss, two-body relaxation, and a realistic treatment of 
the external tidal field. We vary the particle number, the 
initial concentration of the clusters and their orbits to see 
how these parameters affect the dynamical evolution of the 
clusters. Our main focus is to study how quantities which 
can be checked observationally, as e.g. the slope of the mass 
function or the fraction of cluster mass in compact remnants, 
change with time. From our results and observational data, 
we can draw conclusions about the starting conditions of 
individual clusters and the most likely form of the IMF of 
stars in globular clusters. 

The paper is organised as follows: In chapter 2 we de¬ 
scribe our simulation method and the initial set-up of our 
runs. Chapter 3 presents our results for the mass-loss of 
the clusters, the changes in their mass-functions and other 
cluster parameters. In this chapter we will also derive some 
simple fitting formulas which describe the change of these 
parameters as a function of the fraction of time spent until 
complete dissolution. In chapter 4 we will compare our re¬ 
sults with observations of galactic globular clusters and in 
chapter 5 we finally draw our conclusions. 
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2 DESCRIPTION OF THE RUNS 


We simulated the evolution of clusters containing between 
N = 8.192 and 131.028 (128K) stars, increasing the par¬ 
ticle numbers by successive factors of 2. The simulations 
were carried out with the collisional Aarseth N-body code 
NBODY4 (Aarseth 1999) on the recently finished GRAPE6 
boards of Tokyo University. The NBODY4 code uses a Her- 
mite scheme with individual timesteps for the integration 
and treats close encounters between stars by KS (Kustaan- 
heimo & Stiefel 1965) and chain regularisations (Mikkola 
& Aarseth 1990, 1993). NBODY4 was specially adapted to 
make use of the GRAPE6 hardware, which was used to cal¬ 
culate the gravitational forces between the stars. 

Our clusters moved on circular or eccentric orbits 
through an external galaxy which followed a logarithmic po¬ 
tential <j){RG) ~ VqIt^Rg with circular velocity Vg = 220 
km/sec. For the integration, we used a coordinate system 
where the cluster remained at the origin and the center of 
the galaxy moved around the cluster, following the cluster 
orbit. Such a coordinate system represents an accelerated 
but non-rotating one. The force from the galactic center and 
the force due to the cluster acceleration were applied to each 
star when it was advanced on its orbit through the cluster. 

In case of an eccentric orbit, we defined an eccentric¬ 
ity parameter e as e = (Ra — Rp)/{Ra + Rp), where Ra 
and Rp are the apogalactic and perigalactic distances of 
the cluster. Most runs for eccentric orbits were made with 
e = 0.5. Runs were also made with eccentricity parameters 
of e = 0.2, 0.3, 0.7 and 0.8 to study the influence of orbital 
eccentricity on the lifetime of the cluster. 

Our clusters followed King profiles initially with central 
concentrations Wo = 5.0 and 7.0. We also tried Wo = 3.0 
models, but found that they were susceptible to a quick dis¬ 
ruption as a result of the stellar mass-loss. This is in agree¬ 
ment with results of Fukushige & Heggie (1995), who also 
found that King Wo = 3.0 models dissolve quickly unless 
the initial mass-function of the cluster stars is very steep. 
We therefore decided to make no series of runs for Wo = 3.0. 
Cluster radii were adjusted such that the tidal radius of the 
King model was initially equal to the tidal radius of the 
external tidal field, calculated according to 


rt 



R. 


2/3 


( 1 ) 


where rtic is the mass of the cluster, Vg the circular velocity 
of the galaxy and Rq the distance of the cluster from the 
galactic center. In case of an eccentric orbit, the cluster size 
was adjusted such that the tidal radius of the King model 
was equal to the tidal radius at perigalacticon. 

Stars were removed from the simulation if their distance 
from the cluster center exceeded two tidal radii in case of a 
circular orbit, and two apogalactic tidal radii in case of an 
eccentric orbit. All simulations were carried out until com¬ 
plete cluster dissolution. During the runs, the cluster center 
was determined by using the method of Casertano & Hut 
(1985). The number of bound stars and the tidal radius of 
the cluster were then determined iteratively by first assum¬ 
ing that all stars still in the calculation are bound and cal¬ 
culating the tidal radius with eq. 1. In a second step, we 
calculated the mass of all stars inside rt and used it to ob¬ 
tain a new estimate for rt- This method was repeated until 


a stable solution was found. For clusters on eccentric orbits, 
the current distance from the center of the galaxy was used 
to calculate the tidal radius in eq. 1. Lagrangian radii were 
then determined from all stars inside a sphere with radius 
rt around the cluster center. 

The initial mass-function of the cluster stars was given 
by a Kroupa (2001) mass-function with upper and lower 
mass-limits of 15 M© and 0.1 M© respectively. Within our 
mass-limits, this mass-function is a two stage power-law with 
near Salpeter like slope above 0.5 M© and a shallower slope 
below 0.5 M©: 


i{m) dm ' 


m dm , 
m~^'^ dm , 


m < 0.5 M© 
m > 0.5 M© 


( 2 ) 


Kroupa (2001) found this mass-function to be the mean 
mass-function of stars in the galactic disc. The adopted slope 
at the low-mass end is also similar to what Zoccali et al. 
(2000) found for the mass-function of low-mass stars in the 
galactic bulge. Our adopted mass-function leads to an initial 
mean mass of <m>= 0.547 M©. 

Stellar evolution was modeled by the fitting formulae of 
Hurley et al. (2000), which give stellar lifetimes, luminosities 
and radii as a function of initial mass m and metallicity Z 
for all phases of stellar evolution from the zero age main- 
sequence to the remnant stages. To apply the formulae, we 
assumed a metallicity of ^ = 0.001 for the cluster stars. 
Mass lost from stars due to stellar evolution was assumed 
to be immediately lost from the star cluster. During the 
integration, the total mass lost due to stellar evolution was 
summed up since it is an important reference quantity for 
our analysis. 

At the start of the simulation, our clusters did not con¬ 
tain binaries. In our simulations, binaries and higher order 
hierarchies which formed during the integration were treated 
as inert, assuming that no collisions or exchange of mass take 
place and that the stellar evolution of the components is not 
altered by their companions. 

We did not apply kick velocities to stars which became 
neutron stars. Although there is evidence for the existence 
of such kick velocities (Hansen & Phinney 1997), the exact 
form of the kick velocity distribution is still a matter of de¬ 
bate. Applying kick velocities would mean that nearly all 
neutron stars are immediately lost from our clusters since 
the sizes of the kick velocities exceed typical escape veloc¬ 
ities from globular clusters by more than an order of mag¬ 
nitude. This seems to be equally unrealistic, given the fact 
that neutron stars are known to exist in globular clusters. As 
a matter of compromise we have therefore chosen a rather 
low value for the upper mass-limit of our adopted IMF, but 
keep all neutron stars which are formed during the run. 

In order to perform the simulations, clusters were scaled 
from physical units to N-body units, in which the constant 
of gravitation, initial cluster mass and energy are given by 
G = 1, M = 1 and Ec = —0.25 respectively. This scaling is 
straightforward since the total mass of the cluster together 
with the radius of the cluster orbit defines the tidal radius 
and the crossing time in physical units. The stellar evolution 
timescales in N-body units are obtained by requiring that 
their ratio to the crossing timescale must remain constant 
after rescaling the clusters to Y-body units. When compar¬ 
ing our simulations to observations, we assume an age of 
T = 13 Gyrs for the galactic globular cluster system. 
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Table 1. Details of the performed W^-body runs. 


Orbital 

N 

^Sim 

Wo 

Orbit 


Mo 

n 

AMsev 

Roiss 


Tec 


Family 





[pc] 

[Mq] 

[pc] 

[%] 

[MYR] 


[MYR] 

(I) 

8192 

4 

5.0 

Cir 

8500 

4497.1 

24.34 

29.5 

4149 ± 

69 

3329 ± 

50 


16384 

2 

5.0 

Cir 

8500 

9003.3 

30.69 

30.0 

6348 ± 135 

5062 ± 

66 


32768 

1 

5.0 

Cir 

8500 

18407.6 

38.95 

32.1 

9696 


8412 



65536 

1 

5.0 

Cir 

8500 

36223.7 

48.81 

32.1 

15197 


13193 



131072 

1 

5.0 

Cir 

8500 

71236.4 

61.15 

32.1 

23769 


21339 


(ii) 

8192 

4 

7.0 

Cir 

8500 

4401.6 

24.17 

27.8 

3688 ± 

75 

1672 ± 

42 


16384 

2 

7.0 

Cir 

8500 

8927.7 

30.60 

29.8 

5896 ± 

65 

2886 ± 

179 


32768 

1 

7.0 

Cir 

8500 

18013.5 

38.67 

31.1 

9694 


4874 



65536 

1 

7.0 

Cir 

8500 

35611.1 

48.53 

31.9 

15129 


7869 



131072 

1 

7.0 

Cir 

8500 

71699.2 

61.28 

32.9 

25506 


12620 


(HI) 

8192 

4 

5.0 

0.5 

8500 

4480.1 

11.69 

26.5 

1961 ± 

45 

1007 ± 

96 


16384 

2 

5.0 

0.5 

8500 

9021.6 

14.76 

29.0 

3142 ± 

84 

2194 ± 

48 


32768 

1 

5.0 

0.5 

8500 

17966.3 

18.57 

29.9 

4899 


3613 



65536 

1 

5.0 

0.5 

8500 

36156.8 

23.45 

31.2 

7536 


5924 



131072 

1 

5.0 

0.5 

8500 

71385.2 

29.42 

31.5 

11675 


9332 


(IV) 

8192 

4 

5.0 

Cir 

2833 

4441.8 

11.66 

23.5 

1135 ± 

38 

853 ± 

45 


16384 

2 

5.0 

Cir 

2833 

9023.8 

14.76 

26.0 

1892 ± 

70 

1576 ± 

59 


32768 

1 

5.0 

Cir 

2833 

18273.6 

18.68 

28.4 

3120 


2848 



65536 

1 

5.0 

Cir 

2833 

35862.7 

23.39 

28.7 

5092 


4635 



131072 

1 

5.0 

Cir 

2833 

71218.1 

29.40 

29.2 

8324 


7656 


(V) 

8192 

4 

5.0 

Cir 

15000 

4489.1 

35.54 

31.2 

7559 ± 126 

5943 ± 

307 


16384 

2 

5.0 

Cir 

15000 

8808.4 

44.49 

31.1 

11495 ± 412 

8724 ± 

641 


32768 

1 

5.0 

Cir 

15000 

18205.4 

56.67 

33.1 

17262 


14241 



65536 

1 

5.0 

Cir 

15000 

35914.9 

71.02 

32.9 

26868 


22515 



131072 

1 

5.0 

Cir 

15000 

71952.3 

89.59 

33.7 

40212 


36424 


(VI) 

32768 

1 

5.0 

0.2 

8500 

17981.2 

29.48 

30.7 

7834 


6342 



65536 

1 

5.0 

0.2 

8500 

35608.5 

37.04 

31.5 

12383 


9514 


(VII) 

32768 

1 

5.0 

0.3 

8500 

18300.0 

25.73 

31.2 

6709 


5205 



65536 

1 

5.0 

0.3 

8500 

35932.9 

32.22 

32.1 

10762 


8352 


(VIII) 

32768 

1 

5.0 

0.7 

8500 

17957.3 

12.15 

28.8 

3111 


2093 



65536 

1 

5.0 

0.7 

8500 

35749.5 

15.29 

29.8 

4940 


3470 


(IX) 

32768 

1 

5.0 

0.8 

8500 

18025.7 

8.94 

28.1 

2309 


1461 



65536 

1 

5.0 

0.8 

8500 

36116.5 

11.27 

29.6 

3549 


2289 



Notes: 

1. For clusters on eccentric orbits, the value given for Rq is the apogalactic distance 


Table 1 gives an overview of the simulations performed. 
It shows the number of cluster stars N, the number Nsim 
of simulations for each cluster, the initial concentration and 
orbital type, the maximum distance from the galactic center 
and the initial mass and tidal radius of the cluster. Also 
shown are the fraction of mass lost due to stellar evolution, 
and the dissolution and core collapse times of the clusters. 
The dissolution times were defined to be the time when 95% 
of the mass was lost from a cluster. Core collapse times were 
determined from the first minimum of the inner lagrangian 
radii. For models for which more than one run was made. 
Table 1 gives the mean dissolution and core collapse times of 
all runs performed and the standard deviation of individual 
runs around the mean. 

Family (I) represents our standard set of runs which 
are clusters moving on circular orbits with Rq = 8500 pc 
and following density profiles with concentration parame¬ 
ters Wo = 5.0. The other cluster families are variations of 
this which study the influence of the initial condition on the 
cluster evolution. Family (II) contains clusters with higher 
central concentration Wo = 7.0, and the following families 


contain clusters moving in eccentric orbits with Ra ~ Rg- 
The next two families are clusters moving on circular orbits 
at smaller and larger galactocentric distances. The last sets 
of simulations contain clusters on orbits with different eccen¬ 
tricities but otherwise identical initial properties to family 

(I). 

3 RESULTS 
3.1 Dissolution 

We start our discussion by presenting the results for the dis¬ 
solution times. Star clusters dissolve as a result of stellar evo¬ 
lution, two-body relaxation and external tidal shocks, which 
all act on different timescales and scale differently with the 
particle number. Stellar evolution is most efficient directly 
after cluster formation. Its influence on the dynamical evo¬ 
lution of star clusters diminishes thereafter since low-mass 
stars retain a much larger fraction of their initial mass due 
to stellar evolution than high-mass stars (see Weidemann 
2000 and Fig. 18 of Hurley et al. 2000). Two-body relax- 


© 2002 RAS, MNRAS 000, 1-23 



Dynamical evolution of star clusters in tidal fields 5 


8K < N < 128K. = 5.0, R;. = 8.5 Kpc 



T [MYR] 


8K < N < 128K, = 5.0, R;, = 2.833 Kpc 



T [MYR] 


BK < N < 128K, = 5.0, Eccentric orbit 



T [MYR] 


8K < N < 128K, = 7.0, = 8.5 Kpc 



Figure 1. Evolution of the bound mass with time for four different families of models: a) Wq = 5.0, circular orbit with Rq = 8.5 kpc, 
b) Wo = 5.0, eccentric orbit, c) Wq = 5.0, circular orbit with Rq = 2.833 kpc, d) Wq = 7.0, circular orbit. Clusters with 8 K stars are 
shown by a solid line, 16K stars with a dotted line, 32K by a short dashed, 64K by a long dashed, and 128K by a dot dashed line. The 
upper curves in each panel show the stellar mass loss from bound stars, which reaches values around 30% for most cases. 


ation arises from close encounters between cluster stars and 
leads to a slow diffusion of stars over the tidal boundary. It 
acts on a timescale (Spitzer 1987) 


Trh 


r 


■i/2 

h 


<m> G ln( 7 A'') 


(3) 


where th is the half-mass radius of the cluster, N the num¬ 
ber of cluster stars and 7 the Coulomb logarithm. Exter¬ 
nal tidal shocks happen in our simulations during pericenter 
passages of clusters on eccentric orbits. Their strength does 
not change with the particle number and depends only on 
the density profile of the cluster. 
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Fig. 1 depicts the evolution of bound mass with time 
for the first four families of cluster models. For N = 8K and 
16K clusters, the average evolution is shown. The strong de¬ 
crease of bound mass at the start is due to the mass-loss 
of high-mass stars from stellar evolution and levels off af¬ 
ter about 25% of the cluster mass is lost. For clusters with 
Wo = 5.0 there is some indication that the mass-loss from 
stellar evolution is accompanied by some induced mass-loss 
resulting from an overflow of stars over the tidal boundary. 
This mass-loss is much reduced in case of the more concen¬ 
trated Wo = 7.0 models. It can be seen in Fig. 1 and Table 1 
that the total mass lost due to stellar evolution of the bound 
stars is almost the same in all cases, since we usually obtain 
values around 30%. Most discrepant are small-Ai clusters 
moving at galactocentric distances of 77 = 2.8 kpc, whose 
lifetimes are comparable with the timescale of stellar evolu¬ 
tion. If we take a population of stars which initially follow 
our adopted IMF and evolve them with the Hurley et al. 
(2000) stellar evolution routines, we find that 33.1% of the 
mass is lost due to stellar evolution after 2.0 • 10^ Myr. This 
is quite close to the values found in most of our runs. Mass 
loss due to stellar evolution happens therefore very early in 
the evolution and is nearly finished by the time the other 
mass-loss mechanisms become important for the mass-loss. 

While clusters on circular orbits lose mass fairly 
smoothly, clusters on eccentric orbits lose mass in distinct 
steps as stars escape mainly while the cluster is near peri- 
galacticon (see Fig. lb). This is shown in greater detail for 3 
clusters on e = 0.5 orbits in Fig. 2. During each perigalactic 
passage, stars at or beyond the perigalactic tidal radius re¬ 
ceive a strong velocity kick from the varying tidal field and 
most are stripped away from the clusters, causing a decrease 
in bound cluster mass. As the clusters move outward again, 
their tidal radii expand and the clusters can fill a much larger 
volume of space. Clusters will therefore undergo an expan¬ 
sion which is driven by relaxation, although only few stars 
move over the tidal boundary during this phase. The next 
perigalactic passages removes most stars scattered beyond 
the perigalactic tidal radius. One would expect that the frac¬ 
tion of stars removed at each passage should decrease with 
increasing Thh/Tore, or, if the orbital timescale is con¬ 
stant, with the particle number. Fig. 2 shows that this is 
indeed the case. 

Fig. 3 shows the cluster lifetimes as a function of the 
initial mass and orbital type. For the same orbit. King 
Wo = 5.0 clusters have lifetimes very similar to the King 
Wo = 7.0 ones, the difference never exceeding 10%. The 
initial concentration seems to influences the lifetimes only 
very little, provided clusters survive the initial phase when 
mass-loss from stellar evolution dominates. 

The horizontal line in Fig. 3 shows the age of the cluster 
system, which we assume to be 13 Gyrs. At the solar radius, 
clusters with mass less than M — 2.8 • 10^ M© are com¬ 
pletely destroyed within one Hubble time. This value rises 
to 8.1 • 10^ Mq in case of an eccentric orbit with e = 0.5, 
and 1.49 • 10® Mq for a cluster at Rg = 2.833 kpc from 
the galactic center. Several clusters in Table 1 have life¬ 
times well exceeding a Hubble time. The cluster starting 
with N = 128K stars, Wo = 5.0, and moving on a circu¬ 
lar orbit with Rg = 8.5 kpc for example still has a bound 
mass of M = 1.7 • IO^Mq after T — IZ Gyr and would have 
evolved into a small globular cluster. 


0,8 



1000 1500 2000 2500 

T [NBODY] 


Figure 2. Evolution of bound mass for clusters on eccentric or¬ 
bits. Shown is the evolution for clusters with (from bottom to 
top) N = 8K, 16K, 128K stars on e = 0.5 orbits. Solid lines mark 
pericenter passages. During these passages, the number of bound 
stars drops sharply as a result of the shrinking tidal radius. In the 
intermediate phase between two pericenter passages, the number 
of cluster stars changes only slowly. The fractional mass-loss per 
passage is strongest for low-Al clusters. 



Initial cluster mass [M^j 


Figure 3. Dissolution times as a function of the initial cluster 
mass and orbital type. The horizontal dashed line shows an age of 
13 Gyrs. Solid lines show a scaling proportional to , which 

provides a good fit for King Wo = 5.0 clusters. The dotted line 
is a scaling law proportional to Tr^, which fits the lifetimes of 
the King Wo = 7.0 clusters. 
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Figure 4. Lifetimes of clusters moving on orbits with different 
eccentricities e but the same apogalactic distances. Lifetimes are 
divided by the lifetime of a cluster moving on a circular orbit 
with radius equal to the apogalactic radius of the clusters on the 
eccentric orbits. The solid line shows a relation (1 — e), which 
provides a satisfactory fit for all eccentricities. 
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With the exception of clusters moving at Ra = 2.8 kpc, 
we always obtain scaling laws T ~ with x equal to 0.75 
for King Wo = 5.0 clusters moving on similar orbits. King 
Wo = 7.0 clusters show a slightly steeper scaling law with 
X = 0.82. Both values are comparable to what Baumgardt 
( 2001 , 2002 ) found for single-mass clusters in tidal fields, 
and also close to the scaling laws obtained by the ’Collabora¬ 
tive Experiment’ (Heggie et al. 1998) for multi-mass clusters 
without stellar evolution. The dependence of the lifetimes on 
the relaxation time is weaker than linear since potential es¬ 
capers need a considerable amount of time to escape from 
clusters in tidal fields (Fukushige & Heggie 2000). While 
they remain trapped to the cluster, they experience further 
scatterings with other cluster stars, which complicates the 
escape process and leads to a deviation from a linear scaling 
with the relaxation time. 

Following Baumgardt (2001, eq. 16), we expect the dis¬ 
solution time of a cluster to be given by: 

Tdiss = k TTih ^choss • (4) 


Using eq. 3 and the fact that the crossing time is propor¬ 
tional to rf/^ fVG M, we can rewrite this as 


Tdiss = k' 




Gl/2ml/2 In(yiV) / 1^1/2 JVl/2„ll/2 


3/2 


(5) 


where m is the mean mass of a star in the cluster. If we 
assume that all radii in a cluster scale with its tidal radius, 
Th is proportional to Vt, we can rewrite eq. 5 with the help 
of eq. 1 to obtain 


Tdiss = k' 


N 


^ln{'yN) 
or, equivalently 


Rg 


Tdi 


[Myr 


= P 


N 


Rg 


ln{'y N) J [kpc] V 220 km/sec J 


Vg 


y. 


( 6 ) 


(7) 


The exact value of the factor 7 in the Coulomb logarithm 
depends on the mass-spectrum and is rather uncertain for 


Figure 5. Ratio of 7V-body to predicted lifetimes according to 
eq. 10 as a function of the W^-body lifetimes for the same clus¬ 
ters as in Fig. 3. The errorbar gives the uncertainty of a single 
calculation, assumed to be 2%. For most clusters, no systematic 
deviation exists between our prediction and the A'^-body results. 
The exception are clusters dissolving in less than 3 Gyrs, whose 
lifetimes are comparable to the timescale of stellar evolution. 


multi-mass clusters. In addition, it is almost certainly vary¬ 
ing in the simulated clusters since the mass-spectrum it¬ 
self changes during the evolution. Fortunately, due to the 
large number of cluster stars, our results are rather insen¬ 
sitive to the exact value of 7 , so we can adopt the value 
found by Giersz & Heggie (1996) for their multi-mass clus¬ 
ters: 7 = 0 . 02 . From a fit to the lifetimes of the Wo = 5.0 
clusters in circular orbits, we then obtain x = 0.75, /3 = 1.91. 
The corresponding values of x and /? for Wo = 7.0 clusters 
are x = 0.82 and (3 = 1.03. 

Despite the fact that clusters on eccentric orbits are 
subject to external tidal shocks, their lifetimes can also be 
fitted by a scaling proportional to TTiff , similar to the circu¬ 
lar case (see Fig. 3). Hence, in order to predict lifetimes we 
only have to determine how the dissolution time of a cluster 
depends on the eccentricity of its orbit. In order to do so, we 
performed a series of runs of clusters with identical initial 
parameters, but moving on orbits with different eccentrici¬ 
ties e, starting at similar apogalactic distances (families VI 
to IX of Table 1). Fig. 4 depicts how the dissolution times 
of these clusters change as a function of orbital eccentricity. 
The solid line shows a relation according to 

Tdiss{() = Tdiss{Q) ■ 73 ~ P) > ( 8 ) 

which provides a good fit to the V-body results for both 
particle numbers. Here Tdiss{0) is the lifetime of a clus¬ 
ter moving on a circular orbit with radius similar to the 
apogalactic radii of the eccentric orbits. Eq. 8 will be invalid 
when the predicted lifetime of a cluster becomes compara- 
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ble to its orbital timescale. The corresponding relation for 
similar perigalactic distances would be 

Toissie) = Tdiss{ 0 ) • (1 + e) , (9) 


i.e. the lifetime of a cluster on an eccentric orbit can be 
increased by at most a factor of 2 compared to the lifetime of 
a cluster moving on a circular orbit at the same perigalactic 
distance. A similar, relatively small increase was also found 
by Joshi et al. (2001) in Monte Carlo simulations of a cluster 
on an e = 0.6 orbit. 

Summarizing, the lifetimes of clusters starting with a 
Kroupa (2001) IMF and moving on orbits with eccentricity 
e through a logarithmic external potential are given by the 
following formula: 


Tjjiss 


N 


ln{Q.02N) 


" Rg ( Vg 
[kpc] \^220km/sec 


( 1 -£),( 10 ) 


with the values of /? and x depending on the initial concen¬ 
tration. Fig. 5 shows the difference between the A^-body life¬ 
times and our prediction for the clusters of Fig. 3. The typi¬ 
cal uncertainty of a single calculation, based on the values in 
Table 1, seems to be roughly 2%. Most A-body results are 
therefore compatible with our prediction. Most discrepant 
are small-A clusters on Re = 2.8 kpc orbits, which dissolve 
considerably faster than predicted. The reason is that they 
contain high-mass stars for a much larger fraction of their 
lifetime. Since high-mass stars are very effective in scattering 
low-mass stars out of a cluster, the cluster lifetimes become 
smaller. Nevertheless, our prediction is accurate to within a 
few percent for most cases. Chapter 4.1 we will compare our 
formula with some results from the literature and discuss its 
implication for the depletion of the galactic globular cluster 
system. 

We finally discuss the mass-loss rate. Fig. 6 shows how 
the cluster mass changes with time after the mass-loss curves 
were normalized by the mass not lost due to stellar evolu¬ 
tion, i.e. 


Mrel{T) 


M{T) 

M(0) - Msev{T) ’ 


( 11 ) 


where M{T) and Msev{T) are the cluster mass and cumula¬ 
tive mass lost due to stellar evolution at time T respectively. 
Shown is the case Wo = 7.0 which has the smoothest mass- 
loss curve. We have normalized the physical time by the 
lifetime of each cluster calculated according to eq. 7. The 
differences among the different curves are relatively small 
and they are all fairly close to a straight line approximation 
for most part of the evolution. Only near the end, a leveling 
off of the mass-loss can be seen. A linear dependence should 
be a good approximation to the mass-loss curve, since, as 
a cluster loses mass, its tidal radius shrinks proportional to 
rt ~ and, assuming a self-similar behavior, all other 

radii shrink in a similar way. Since the relaxation time de¬ 
pends on r and M roughly as Trh ~ \/M Trh is 

proportional to the mass left in the cluster. From the results 
obtained so far, we would expect the mass-loss rate to be 
given by dM/dT ~ —M/Tre = which varies only 

slowly with M since x is not too far from unity. 

Hence, if we assume that the mass-loss due to stellar 
evolution takes place immediately after cluster formation, 
and that later on globular clusters lose mass linearly, the 



T/TDISS 


Figure 6. Evolution of mass with time for the Wq = 7.0 clusters 
after the mass lost due to stellar evolution has been subtracted 
from the total mass loss and the timescales have been divided 
by the dissolution times from eq. 7. Clusters lose mass roughly 
linearly until shortly before the end. 


mass still bound to a globular cluster with initial mass Mq 
can be approximated as 

M(T) = 0.70 Mo (1 - T/Toi^s) , (12) 

for T < Toiss- 


3.2 Evolution of the mass-function 

Fig. 7 depicts the evolution of the combined mass-function 
of main-sequence stars and giants as a function of time. 
The bend in the Kroupa (2001) IMF at m = 0.5 M© can 
clearly be seen. It survives until fairly late but is difficult 
to distinguish by the time 90% of the mass is lost. Below 
m = 0.5 Mq, the mass function can be characterised by 
a single power-law throughout most part of the evolution. 
Only for clusters on eccentric orbits, some deviations from 
a power-law occur near the end of the lifetime at the low- 
mass end. There are indications that the mass-function of 
some galactic globular clusters cannot be characterised by 
single power-law mass-functions (de Marchi, Paresce & Pu- 
lone 2000). If true, these features are likely to be primordial, 
since Fig. 7 shows that dynamical evolution does not lead 
to the formation of such features. 

Fig. 7 also shows the preferential depletion of low-mass 
stars as a result of mass-segregation and overflow over the 
tidal boundary. The depletion already sets in before 30% 
of the lifetime is over. After about 60% of the lifetime has 
passed, low-mass stars are depleted to such an extent that 
the number of cluster stars starts to drop towards the low- 
mass end. 

Differences between large and small-A models are rela¬ 
tively small for any orbital family if we compare the mass- 
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Figure 7. Combined mass function of main-sequence stars and giants for different orbital families. Shown are the mass-functions at 
the start, and when 30%, 60% and 90% of the cluster lifetime are over. Solid lines are for the biggest simulated clusters of each family 
{N = 128K), dashed lines show the sum of the four smallest clusters with N = 8K stars. Small-clusters were shifted to contain the 
same number of low mass stars (m < 0.5 Mq) than the high-A'^ runs. The preferential depletion of low-mass stars leads to mass-functions 
which decrease towards small m towards the end. The details of this depletion are nearly independent of the initial particle number and 
orbital type. 


functions at the same fractional lifetime, despite the fact 
that there is already a factor of 16 difference in the parti¬ 
cle numbers. It should therefore be safe to apply our results 
also to globular clusters, which mostly have particle num¬ 
bers within a factor of 10 of our largest runs. Furthermore, 


differences between clusters of different orbital families are 
also small, indicating that the mass-function changes in an 
universal way, independent of particle number and orbital 
type of the parent cluster. It is often stated in the litera¬ 
ture that the preferential depletion of low-mass stars is due 
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Figure 8. Slope of the mass function of low-mass stars as a func¬ 
tion of the fraction of time spent until complete cluster dissolu¬ 
tion. a) Slopes for individual W^-body runs, b) Differences between 
the individual slopes and an average taken from all runs shown 
in panel a). The differences are small and could entirely be due 
to finite-Af scatter, c) Difference of the mean from a least-square 
fitting curve according to eq. 13. 


to mass-segregation and external tidal shocks. However, our 
simulations show that relaxation is responsible for the de¬ 
pletion of low-mass stars, since we get the same depletion 
whether clusters are subject to external tidal shocks (as in 
Fig. 7b) or not. 

Fig. 8 depicts the evolution of the slope of the mass 
function at the low-mass end in dependence of the time 
elapsed until cluster dissolution. The slope was determined 
from a fit to the distribution of stars with masses m < 
0.5 Mq. Note that we have put stars into bins equally spaced 
in log m, so the slope x of the mass-function is by one smaller 
than in case of a plot over m, i.e. x = a — 1. Fig. 8 confirms 
that the slope of the mass-function changes in all clusters 
in more or less the same way. During the first 20% of the 
cluster lifetime, the mass-function remains nearly constant, 
which might arise from the fact that in the beginning masses 
of stars with m < 0.5 Mq are fairly small compared to the 
mass of the most heavy stars in the cluster, hence they will 
be depleted at similar rates. Later, the slope changes with 
increasing speed since low-mass stars become more and more 
depleted. Differences between individual runs (Fig. 8b) re¬ 
main small and could entirely be due to statistical errors, 
which are of order Aa; = 0.02 at the start of the runs and 
increase as the number of cluster stars decreases. 

We can fit the mean evolution by a formula like: 


X = 0.3-1.51 



2 

-bl.69 



-1.50 



where we have chosen the functional form first and deter¬ 


mined the coefficients from a least mean-square fit to the N- 


Figure 9. Same as Fig. 8 but now in dependence of the mass 
fraction which is still bound to the clusters. 

body results. The differences between eq. 13 and the mean 
curve from the A-body results are small and are shown in 
Fig. 8c. Eq. 13 should therefore give an accurate prediction 
for the change of the stellar mass-function in a globular clus¬ 
ter, provided its dissolution time is known. 

Similar results can also be found for the dependence of 
the slope of the mass-function on the fraction of mass that 
is already lost from a cluster (Fig. 9). In this case, the slope 
does not change significantly until nearly half the cluster 
mass is lost since in the beginning most mass is lost due to 
stellar evolution. We can fit the mean evolution by: 

a: = 0.3 - exp(0.67 - 6.19M/Mo - 3.24 (M/Mof) (14) 

The lower panels of Fig. 9 show the difference between the 
mean curve and individual runs and the difference between 
a mean curve and our fitting function. Differences are of the 
same size as before. 

Mass-segregation changes the local mass-function rel¬ 
ative to the global one since stars heavier than the mean 
sink towards the center and fighter ones move outward. If 
uncorrected for, this could bias the observed mass-function, 
as observations usually cover only a limited range in radii, 
most often around the half-mass radius. In order to deter¬ 
mine to which extent mass-segregation changes local mass- 
functions, we calculated local mass-functions from all stars 
with m < O.5M0 at different radial shells and compared 
them with the global one. Fig. 10 shows the results for the 
N — 128K star cluster from our standard family. Shown are 
differences for 3D shells (panel a) and (projected) 2D shells 
(panel b). In order to allow the most easy comparison with 
observations, projected shells consider the emitted fraction 
of cluster light. 

In the beginning, the mass-functions are the same at 
any radius since our clusters start without mass-segregation. 
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Figure 10. Difference of the local mass-function minus the global one at four different times for the run with N = 128K stars and 
Rq = 8.5 kpc. Shown are differences at the start of the calculation (solid line) and when 30% (dotted), 60% (short dashed) and 90% (long 
dashed) of the lifetime have passed. The left panel shows differences for 3D mass-shells, the right panel shows differences in projection 
and as a function of the enclosed fraction of cluster light. Inside the half-mass radius, local mass-functions are steeper than the global 
one as a result of mass-segregation. Throughout the run, mass-segregation is increasing in the inner parts. Differences in projection are 
smaller than the 3D ones. Best agreement with the global mass-function is achieved just outside the half-light radius. 


Mass-segregation is however quickly established and one 
finds more negative slopes for the mass-function in the in¬ 
ner parts since heavier stars sink into the center and lighter 
stars move outwards. The largest differences are found in the 
center and amount to about Ax = 1.0 in the 3D case close 
to the end. The reason for the strong difference compared 
to the average slope over all radii is that low-mass stars are 
almost completely absent in the cluster center. Outside the 
half-mass radius, local mass-functions differ much less from 
the global one and are always shallower. As one would ex¬ 
pect, differences are smaller if viewed in projection (panel 
b). The radii containing between 50-80% of the emitted clus¬ 
ter light show the smallest difference from the global mass- 
function. For the cluster shown, this amounts typically to 
distances of 1 to 1.8 (projected) half-light radii. Similar re¬ 
sults are obtained for the other clusters. Mass-segregation 
is stronger and sets in faster in clusters undergoing core¬ 
collapse earlier in their evolution, as e.g. clusters starting 
from Wo = 7.0 models. Here differences between local and 
global mass-functions are about Ax = 0.2 larger in the cen¬ 
ter than in the Wo = 5.0 case. 


3.3 Cluster composition and mass-to-light ratio 

The question to which extent globular clusters contain dark 
matter is still a matter of debate (see Heggie & Hut 1996). 
Dark matter could be present either in the cluster centers 
in the form of intermediate mass black holes (Gerssen et al. 
2002, Gebhardt, Rich & Ho 2002), or in the form of faint 
stellar mass white dwarfs and neutron stars. Previous studies 


(Vesperini & Heggie 1997, Giersz 2001) have shown that 
the fraction of white dwarfs is rising in a star cluster as it 
evolves. According to Giersz (2001), more than half the total 
cluster mass can be in the form of white dwarfs close to the 
end of its lifetime, depending on the initial mass-function 
and cluster lifetime. Since our initial mass-function and the 
results for the cluster lifetimes are different from previous 
studies, it might be interesting to look how the mass fraction 
of compact remnants changes in our clusters. 

Fig. 11 depicts for two different runs the fraction of 
mass in different stellar groups as a function of the number 
of stars left in the cluster. We have divided main-sequence 
stars and giants into two categories, high and low-mass stars, 
depending on whether their mass exceeds 0.7 M© or not. In 
the beginning, most mass in the cluster can be found in 
high-mass main-sequence stars. Their mass-fraction drops 
steadily since stellar evolution constantly transforms the 
high-mass end of the mass-function into compact remnants. 

The mass-fraction of low-mass stars increases initially 
due to the mass-loss of the high mass stars. Later, their 
mass-fraction decreases since they are the lightest mass com¬ 
ponent and are therefore preferentially removed from the 
cluster. 

The decrease in the mass-fraction of main-sequence 
stars is compensated for by a constant increase of the mass- 
fraction of white dwarfs and, very close to the end, also by an 
increase in the mass-fraction of neutron stars. For the clus¬ 
ter with N = 128K stars initially, white dwarfs become the 
dominant mass group after roughly 80% of the cluster stars 
are lost. The mass and number fractions of giants are always 
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Figure 11. Fraction of mass in different stellar groups in depen¬ 
dence of the number of stars still bound to the cluster. Solid lines 
show the evolution for the run with N = 128K stars, Wq = 5.0 
and Rc = 8.5 kpc. The dashed lines show the evolution for a 
cluster with N = 32K stars but otherwise similar initial condi¬ 
tion. The fraction of compact remnants increases throughout the 
evolution. The mass-fraction of low-mass stars increases in the 
beginning due to the mass-loss of high-mass stars. Later their 
fraction decreases as a result of the preferential depletion of low- 
mass stars. 


very small and they are therefore not shown separately. The 
dotted lines show the evolution for a cluster with a smaller 
initial particle number. It can be seen that the overall evo¬ 
lution is similar, except that the fraction of white dwarfs 
is smaller since in a cluster with shorter dissolution time a 
smaller fraction of stars is turned into white dwarfs due to 
stellar evolution. 

For a quantitative understanding of how the fraction 
of mass belonging to the different mass-groups changes, one 
hrst has to know the rate with which stellar evolution turns 
mass into compact remnants. Fig. 12 a) depicts how the 
combined mass-fraction of white dwarfs and neutron stars 
increases with time if we take a population of 500000 stars 
that follow a Kroupa IMF initially and evolve them with the 
Hurley et al. (2000) stellar evolution routines. For T > 10 
Myr, we can fit the whole evolution rather accurately by the 
following formula: 

FcrIsev = 0.024 - 0.048 a -b 0.0241 - 0.00055 2^ 

-b0.00033 2"^ , (15) 

where 2 = log]^g(r). Here T is measured in Myrs and the 
coefficients were determined by a least square ht to the evo¬ 
lution of the test population. The differences between our 
fitting formula and the real evolution always remain smaller 
than AF = 0.005. 

Fig. 13 shows the mass-fraction of compact remnants 
for five individual cases, drawn from runs with both small 




Figure 12. Evolution of the mass-fraction of compact remnants 
(left) and the cluster mass-to-light ratio (right) for a population of 
500000 stars following a Kroupa (2001) IMF and evolving under 
the influence of stellar evolution alone. Solid lines represent results 
of simulations, dashed lines the fit according to eqs. 15 and 17. 



Figure 13. Fraction of mass in compact remnants for 4 individual 
runs. From top to bottom: N = 128K, Wq = 5.0, Rq = 8.5 kpc; 
N = 64K, Wo = 7.0; N = 128K, eccentric orbit; N = 32K, Wo = 
5.0, Rg = 2.8 kpc. Solid lines show the evolution before, dashed 
lines after the mass-fraction due to stellar evolution (eq. 15) is 
subtracted from the runs. After this subtraction, the evolution is 
rather similar in all cases. 


and large numbers of stars. There are considerable differ¬ 
ences in the mass-fractions of individual runs arising from 
differences in the dissolution times: Small-Af clusters close to 
the galactic center dissolve considerably faster (in Myrs) and 
turn a much lower fraction of bound stars into white dwarfs 
than clusters at larger Rg- Except for the small-Af cluster 
close to the galactic center, we always And that the mass 
fraction of compact remnants is near 90 % close to the end 
of the cluster lifetime. After we subtract the increase due 
to stellar evolution according to eq. 15 from the different 
cases, the differences become rather small (dashed lines). In 


© 2002 RAS, MNRAS 000, 1-23 








Dynamical evolution of star clusters in tidal fields 13 


agreement with Vesperini & Heggie (1997), we find that the 
mass-fraction of compact remnants is still increasing since 
low-mass main-sequence stars are preferentially lost from 
the clusters. From Fig. 13, we obtain the following fit to the 
overall evolution: 

For = Fbfllggy-1-0.18 (T/Tdiss)^ + 0.16 {T/Td issf'-{16) 

A similar analysis can be done for the cluster mass-to- 
light ratios. Fig. 12 b) depicts the evolution of the M/L 
ratio for a population of stars which evolves due to stellar 
evolution alone. Clusters start with an M/L ratio of 0.08. 
This value increases constantly as the most luminous stars 
disappear. After 13 Gyrs, the M/L ratio has reached 2.0, 
which is within the range of mass-to-light ratios observed 
for galactic globular clusters (Mandushev et al. 1991, Pryor 
& Meylan 1993). For T > 100 Myr, we can fit the evolution 
of the MLR by the following equation: 

M/L\gi^y = 3.632-5.6282-b3.3642^ -0.9192®-b0.l2'‘ ,(17) 

where 2 is again 2 = logjg(r). The scatter around the mean 
relation is somewhat larger than before since the cluster lu¬ 
minosity is dominated by relatively few giants. 

Fig. 14 a) shows the evolution of the M/L ratios for four 
different clusters. In general, all mass-to-light ratios are in¬ 
creasing, although with considerable differences between in¬ 
dividual runs. Differences arise since clusters which dissolve 
quicker do so at an age in which they still contain many lu¬ 
minous stars. Panel b) of Fig. 14 shows the M/L ratios for 
the same clusters after the effect of stellar evolution due to 
eq. 17 has been subtracted: 

AM/L = M/L - M/L\g^y . (18) 

The remaining change is due to the dynamical evolution 
alone. It can be seen that dynamical evolution leads to a 
decrease of the M/L ratios because faint, low-mass stars are 
preferentially lost from the clusters. On the other hand the 
fraction of compact remnants which contribute almost noth¬ 
ing to the cluster luminosity is increasing, so that close to 
the end the M/L ratios increase again. For galactic globular 
clusters, we would predict that clusters in an advanced stage 
of their evolution have smaller mass-to-light ratios than clus¬ 
ters of similar age which have lost only small part of their 
mass, the maximum differences being about AM/L — 0.5. 
Mandushev et al. (1991) and Pryor & Meylan (1993) have 
indeed found a correlation of the average M/L ratio of a 
cluster with its mass in the sense that low-mass clusters 
have smaller M/L ratios. This could be the result of the 
dynamical evolution since, for example, any cluster within 
a few kpc from the galactic center which has a current mass 
M < 10® M 0 must have lost most of its initial mass and 
should be close to final dissolution (see Fig. 25). 

3.4 Mean masses and central remnant fraction 

Fig. 15 depicts the evolution of the mean masses of all stars 
in our standard run. In the beginning, mean masses drop 
throughout the cluster due to the mass-loss of the most 
massive stars from stellar evolution. For the inner radii, 
this is soon compensated by mass segregation, which causes 
heavy-mass remnants and main-sequence stars to sink to¬ 
wards the cluster center. As a result, the mean mass in 
the core rises until core-collapse, when it reaches a value 



Figure 15. Evolution of the mean mass of all stars with time 
for the run with N = 128K stars. Wo = 5.0 and Rq = 8.5 
kpc. Shown are the mean mass of all stars in the cluster (dashed 
line) and (from top to bottom) mean masses in lagrangian shells 
containing 1%, 5-10%, 40-50%, 70-80% and 90-100% of the bound 
mass. The time of core collapse is marked by an arrow. The mean 
mass of stars near the half-mass radius is always close to the mean 
mass of all stars in the cluster. 


of <m>= 1.2 Mq. This is similar to the masses of massive 
white dwarfs and neutron stars and shows that in the cen¬ 
ter compact remnants must make up a significant fraction 
of the stellar population by this time. After core-collapse, 
the mean mass in the center remains nearly constant. Out¬ 
side the half-mass radius, mean masses decrease initially as 
a result of stellar evolution and the drift of massive stars 
into the core. Later, they increase since low-mass stars are 
preferentially lost from the cluster. In the end, mean masses 
are almost the same throughout the cluster. The mean mass 
around the half-mass radius stays close to the mean mass of 
all stars in the cluster (dashed line). 

Fig. 16 shows how the mass-fraction of white dwarfs and 
neutron stars changes in the core as a function of the number 
of stars left in the cluster. Here, the core was defined as the 
cluster part which contains the innermost 2 % of the mass. 
The mass-fraction of remnants increases faster in the center 
than in the rest of the cluster since white dwarfs and neu¬ 
tron stars sink towards the center due to mass-segregation. 
By the time of core-collapse, the central remnant fraction 
has reached 95%, despite the fact that their overall mass- 
fraction is only 60% (see Fig. 11). At core-collapse 40% of 
the central mass is in the form of neutron stars. This value 
stays nearly constant in post-collapse, since most neutron 
stars have masses of order 1.2 M©, similar to the mass of 
the most massive white dwarfs. Hence, the relative mass- 
fraction of both components should not be altered much by 
the dynamical evolution of the cluster. 

Most other clusters also show large central remnant 
fractions of order 90% or more. The exact value depends 
on the time at which core-collapse happens: Clusters which 
undergo core-collapse later and have turned a larger frac¬ 
tion of their stars into compact remnants show generally 
also larger central remnant fractions. 

Lugger et al. (1995) have measured U-band CCD sur- 
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Figure 14. Cluster mass-to-light ratios for the same clusters as in Fig. 13. Panel a) shows the full evolution, panel b) shows the remaining 
part after the effect of stellar evolution has been subtracted. Dynamical evolution decreases the average M/L ratio by about 0.3 to 0.7 
at T/Tbiss = 0-8 and increases it towards the end of the lifetime. 


face brightness profiles for 15 globular clusters out of the 
sample of 21 clusters identified to be core collapse clusters 
by Djorgovski & King (1986). They found that 13 of them 
have power-law slopes in their centers with a mean slope of 
a(r) = r-°“i. 

Fig. 17 shows the mass density profile of the N = 128K 
star, e = 0.5 run at core-collapse time. We have chosen this 
cluster since it has a large number of stars remaining at 
core-collapse and also a core-collapse time which is close to 
the present time. It might therefore be one of our best ana¬ 
logues to a large galactic globular cluster. The strong central 
concentration of white dwarfs and neutron stars can clearly 
be seen. Inside 0.1 pc, their mass-density is more than a 
factor of 30 higher than that of main-sequence stars. The 
density profile of compact remnants can be approximated 
by a power-law with slope p ~ from the centre out to 

about 2 pc. This slope is in good agreement to what Taka- 
hashi et al. (1997) found from Fokker-Planck simulations of 
multi-mass clusters without mass-loss. The main-sequence 
stars follow a shallower profile then the compact remnants 
inside 0.5 pc with p ~ Outside this radius, their den¬ 

sity profile decreases is even shallower up to 3 pc. Beyond 
this radius the profile steepens again due to the presence of 
the tidal field. Since the central density of main-sequence 
stars is much smaller than that of the compact remnants, 
core-collapse is done mainly by the remnants and not by the 
main-sequence stars. Projected and number density distri¬ 
butions show a similar behavior. For the projected number 
densities we obtain slopes of ctcr ~ r~^'^ and ctms ~ 
for compact remnants and main-sequence stars respectively. 
Our density profile for the main-sequence stars is therefore 
in very good agreement to what Lugger et al. (1995) found 
for the mean projected density distribution. Furthermore, 


several clusters in their paper, e.g. NGC 6284, 6325, 6397 
and 6522, also show a saddle like profile at intermediate radii 
between the core and the halo. This could be an indication 
that in these clusters too most mass in the centre is in the 
form of white dwarfs and neutron stars. In contrast, main- 
sequence stars in the run with Wo = 7.0 andA^ = 128K 
stars follow a power-law profile after core-collapse from the 
center all the way up to several pc. A discussion where the 
properties of this cluster are compared with the center of 
the globular cluster M15 and the implications for a possible 
detection of a intermediate mass black hole in this cluster 
(Gerssen. et al. 2002) can be found in Baumgardt et al. 
( 2002 ). 

A close inspection shows that most white dwarfs in the 
core tend to be old, high-mass white dwarfs which are al¬ 
ready very faint. In addition, the neutron stars, if not mem¬ 
bers in close binary systems, would also be difficult to detect. 
In principle, a centrally concentrated distribution of faint 
stars could appear as an intermediate mass central black 
hole to observers, since they also give rise to an increasing 
velocity dispersion of the main-sequence stars. In order to 
test if the compact remnants in the core could appear as a 
black hole, we had a look at the central mass-to-light ratio, 
since a large central mass-to-light ratio could be taken as 
an indication for the presence of a massive black hole in the 
cluster center. 

Fig. 18 shows the evolution of the central mass-to-light 
ratio for the N — 128K star, e = 0.5 run. Here we dehned as 
core the radius of the shell which contains the innermost 2% 
of main-sequence stars in projection. Since high-mass stars 
sink preferentially into the core, the central mass-to-light 
ratio is in the beginning smaller in the core than in the rest of 
the cluster. This changes as more and more white dwarfs and 
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Figure 16. Evolution of the central mass fraction of compact 
remnants for the same run as in Fig. 15. Since mass-segregation 
causes compact remnants to sink into the center, their mass frac¬ 
tion is increasing faster in the center than in the rest of the clus¬ 
ter. By the time of core-collapse (dotted line), the center is almost 
entirely made up out of neutron stars and white dwarfs. 



r [pc] 


Figure 17. Mass density at core collapse for the run with e = 0.5, 
N = 128K stars. The combined mass density of white dwarfs and 
neutron stars is shown by crosses, the density profile of main- 
sequence stars by dots. The central density of compact remnants 
exceeds that of main-sequence stars by more than a factor of 30 
in the core. 



Figure 18. Total and central mass-to-light ratios for the same 
run as in Fig. 17. Central mass-to-light ratios were calculated for 
the shell containing 2% of all main-sequence stars in projection. 
The arrow marks the time of core-collapse. After core-collapse, 
central mass-to-light ratios can reach values of 10 or more due to 
the high fraction of compact remnants in the core. 

neutron stars accumulate in the core, and after core-collapse 
the central mass-to-light ratio reaches values of 10 or more. 
This could explain the high central M/L ratios found for 
some globular clusters: D’Amico et al. (2002) for example 
found that an M/L ratio of 10 is necessary to explain the 
pulsar timings in NGC 6752. Since this cluster has a light- 
profile similar to a post-collapse cluster (Djorgovski 1993), 
the large mass-to-light ratio could be due to a concentration 
of compact remnants in the cluster core. 

A major uncertainty of the present analysis is that our 
clusters did not contain primordial binaries. Since at least 
some fraction of these binaries would contain main-sequence 
stars with masses of m « 0.7Mq, the total mass of the bi¬ 
nary would be comparable to, or even exceed the masses of 
high-mass white dwarfs, so these binaries could also sink to¬ 
wards the core, thereby decreasing the central mass-to-light 
ratio. On the other hand, core-collapse can be achieved only 
if the binary fraction is decreased considerably from the ini¬ 
tial value so that there is no significant energy release from 
binaries any more. In addition, any main-sequence star in 
a binary which has reached the center is likely to be re¬ 
placed by a more massive compact remnant through dy¬ 
namical interactions, and would subsequently be expelled 
from the core. Central mass-to-light ratios might therefore 
be similar to the present values even if clusters started with 
a significant fraction of primordial binaries. 

3.5 Energy equipartition 

Any self-gravitating system which contains stars of different 
masses will show a tendency towards energy partition, such 
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Figure 19. Time evolution of the energy equipartition parame¬ 
ters P for 5 different lagrangian shells for the run with N = 128K 
stars, Wq = 5.0 and Rq = 8.5 kpc. Shown are parameters for 
the shells containing (from bottom to top) 5%, 10-20%, 30-40% 
50-60% and 80-90% of the mass. Energy equipartition is reached 
at the center at core-collapse time (marked by an arrow). Outside 
the half-mass radius, it is not reached within the lifetime of the 
cluster. 


that at any radius heavy stars have the same kinetic energies 
as lighter ones. If stars start with the same average velocities 
initially, this implies that heavy stars sink towards the center 
and lighter stars will move outwards. If the density of high- 
mass stars in the cluster center becomes too large, they will 
form a self-gravitating system of their own, in which case 
equipartition of energies cannot be reached (Spitzer 1969). 

In order to test for equipartition of energies in our runs, 
we define a correlation parameter between the energy of a 
star and its mass as 

p_ T,i (<'!'>-B) {<'^oe rn>-log mi) 

\/Ei(<^> \/Ei(< log"* > - log mi)2 

where the sums run over all stars within a radial shell and 
< T > and < log m > are the mean kinetic energy and log¬ 
arithm of the mass of stars in the shell. The logarithm of 
the mass was chosen in order to give less weight to the few 
stars at the high-mass end. P will be zero if kinetic energy is 
uncorrelated with mass, i.e. energy equipartition is reached. 

Figs. 19 and 20 show the evolution of the P parameters 
for 5 different lagrangian shells for runs with initial central 
concentrations of ITo = 5.0 and ITo = 7.0 respectively. The 
arrows mark the time of the first core-collapse in both runs. 
Energy equipartition is reached for the inner shells at around 
core-collapse time. For radii near the half-mass radius, en¬ 
ergy equipartition is reached only close to the end of the 
lifetime and for the halo it is never reached. This is the case 
even in the run with Wo = 7.0, despite the fact that this run 
spends nearly half its lifetime in the post-collapse phase. We 
find this result for all our runs and similar results were also 
obtained by Giersz & Heggie (1996) in their simulations of 
small-clusters without stellar evolution. We conclude that 
for globular clusters, energy equipartition, if not present in 
some way already initially, can be expected to exist only in 


Figure 20. Same as Fig. 19 but now for the cluster with initial 
concentration Wq = 7.0 and N = 128K stars. 


the centers of clusters which went through a core collapse 
phase. 


3.6 Velocity anisotropy 


It has long been known that isolated star clusters become 
radially anisotropic after core-collapse due to the scatter¬ 
ing of stars out of the cluster core (Spitzer 1987). In a tidal 
field, the situation is complicated by the outflow of stars 
over the tidal boundary. Giersz & Heggie (1997) found in 
V-body simulations that star clusters in tidal fields remain 
isotropic during their evolution except for the outer parts 
which become tangentially anisotropic due to the preferen¬ 
tial loss of stars on radial orbits. Similar results were also 
found by Takahashi, Lee & Inagaki (1997) and Takahashi & 
Lee (2000) in their anisotropic Fokker-Planck simulations. 
There is also mounting observational evidence that at least 
some globular clusters, like e.g. lo Gen (van Leeuwen et al. 
2000 ), have anisotropic velocity distributions, although in 
this case it is not clear whether such velocity distributions 
were primordial, or formed as a result of the dynamical evo¬ 
lution. 

Since it is difficult to incorporate all effects of an ex¬ 
ternal tidal field into a Fokker-Planck code and since most 
V-body simulations performed contain only relatively few 
stars, it might be interesting to look for signs of velocity 
anisotropies in our runs too. Following Binney & Tremaine 
(1986), we define an anisotropy parameter /3 as 


13=1- 


Vt 

2 Vr 


( 20 ) 


where Vr and Vt are the radial and tangential velocity dis¬ 
persions. They are measured in a non-rotating coordinate 
system in which the cluster center is always at rest at the 
origin. Fig. 21 shows the evolution of the velocity anisotropy 
for the run with N = 128K stars from our standard family. 
It can be seen that the outer parts become rapidly tangen¬ 
tially anisotropic. Inspection of the orbits shows that most 
stars in the outer parts are counterrotating, probably due 
to the fact that retrograde orbits are more stable against 
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Figure 21. Evolution of the velocity anisotropy with time for 
different radii for the run with N = 128K stars, Wo = 5.0 and 
Rq = 8.5 kpc. Arrows mark the core-collapse time. Tangential 
anisotropy develops rapidly in the outer parts. Its size remains 
more or less constant until half-mass time, after which it decreases 
slowly. In the inner parts, clusters stay isotropic throughout the 
evolution. 

escape than prograde orbits (Keenan & Innanen 1975). As a 
consequence, the clusters also start to rotate in their outer 
parts. Fig. 21 shows that the tangentially anisotropic veloc¬ 
ity dispersion is restricted to the outermost radii. We find 
a similar behavior for the counterrotation of the stars. We 
note that the amount of tangential velocity anisotropy and 
rotation in the outer parts would have been larger if, as is 
frequently done, we would have chosen a rotating coordinate 
system as reference in which the direction from the cluster 
to the galactic center remains fixed. 

The amount of tangential anisotropy remains more or 
less constant until near core-collapse, when it starts to de¬ 
crease. An explanation could be that the relaxation time of 
a cluster decreases as it loses mass and that a smaller relax¬ 
ation time in the halo implies a more efficient isotropization 
of the stellar orbits. In addition, more stars are scattered out 
of the core into the outer cluster areas in the post-collapse 
phase. Since these stars will be on near radial orbits initially, 
they will also decrease the amount of tangential anisotropy 
in the halo. For clusters in eccentric orbits we obtain simi¬ 
lar results, except that the amount of tangential anisotropy 
is smaller compared to the circular case (Fig. 22) Globu¬ 
lar cluster might therefore also have slightly tangentially 
anisotropic velocity distributions in their halos. 

Except for the outermost radii, clusters are nearly 
isotropic and there is no sign of the radially anisotropic 
velocity dispersion which forms in isolated clusters. These 
results confirm earlier results obtained by Takahashi & Lee 
(2000) and Giersz & Heggie (1997). 

Takahashi & Lee (2000) found that there is a difference 


Figure 22. Same as Fig. 21 for a cluster in an eccentric orbit 
with e = 0.5. Tangentially anisotropic velocity distributions still 
form in the outer parts, but the amount of anisotropy is reduced. 

in the mean anisotropies of high and low-mass stars, espe¬ 
cially in the cluster halos. In their simulations, high-mass 
stars always followed more radially anisotropic velocity dis¬ 
tributions than low-mass stars, the difference depending on 
the fraction of mass lost from a cluster and the distance from 
the cluster center. Fig. 23 shows the anisotropy parameter (3 
for 4 different groups of main-sequence stars at several times 
during the cluster evolution. Close to the end, the velocity 
distributions become indeed different. In the halo we find 
differences of up to = 0.3, while near the cluster center 
differences remain small at any time. The differences could 
be due to mass-segregation: Since massive stars are more 
centrally concentrated, they tend to be located at lower en¬ 
ergy levels which have more isotropic velocity distributions. 
Since mass-segregation itself develops only fairly late in the 
halo, it also takes very long until there are differences in the 
anisotropies. 


4 COMPARISON WITH GALACTIC 
GLOBULAR CLUSTERS 

4.1 Dissolution times and the depletion of the 
galactic globular cluster system 

Dinescu et al. (1999) have compiled a catalog of absolute 
proper motions for 38 globular clusters. Besides the cluster 
orbits, they also determined destruction rates for individ¬ 
ual clusters, following the formalism of Gnedin & Ostriker 
(1997). Their lifetimes provide a good check for our formula 
for the lifetime. In order to compare the lifetimes, we first 
took the orbital parameters for each cluster from Table 5 of 
Dinescu et al. (1999) and calculated the present day cluster 
masses Me from the absolute magnitudes in the McMas- 
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Figure 24. a) Lifetimes from Dinescu et al. (1999) against lifetimes derived in this work for a sample of 38 globular clusters. There 
is good agreement except for clusters close to dissolution for which we predict slightly longer lifetimes, b) Ratio of the lifetimes as a 
function of the relative importance of tidal shocking compared to the total destruction rate as given by Dinescu et al. (1999). For clusters 
for which shocking is important, Dinescu et al. (1999) derive smaller lifetimes. 


ter database of globular cluster parameters (Harris 1996), 
assuming a mass-to-light ratio of M/Lv ~ 2. We then cal¬ 
culated iteratively the initial mass Mo necessary to produce 
a cluster with mass Me after a Hubble time and the remain¬ 
ing lifetime for each cluster, using eqs. 10 and 12. 

Fig. 24 compares our results for the remaining lifetime 
with the destruction rates obtained by Dinescu et al. (1999) 
for the galactic model of Paczyriski (1990). There is quite 
good agreement between the lifetimes determined from N- 
body simulations and those determined from Fokker-Planck 
calculations. The exception might be clusters close to de¬ 
struction, where the lifetimes of Dinescu et al. (1999) seem 
to be somewhat smaller than ours. Panel b) compares the 
differences in the lifetimes as a function of the importance of 
disc and bulge shocking for the total destruction rate. Both 
parameters are taken from Dinescu et al. (1999). It can first 
be seen that the clusters for which tidal shocks are least im¬ 
portant have longer lifetimes according to our formula. Close 
inspection of the orbits shows that nearly all of these clusters 
are orbiting at galactocentric distances larger than 10 kpc, so 
the differences could be due to differences in the underlying 
galactic models. Secondly, clusters which dissolve mostly due 
to external tidal shocks also have longer lifetimes according 
to our formula. This could be a hint that our neglect of disc- 
shocking overpredicts the lifetimes in these cases. However, 
since most of these clusters are also close to dissolution, it is 
difficult to determine what is responsible for the difference. 
Nevertheless, the overall agreement is reasonably good and 
we find the same agreement for the lifetimes obtained by 
Dinescu et al. for a different galactic model. 

We finally note that the lifetimes of eq. 10 are also in fair 
agreement with those obtained by Vesperini & Heggie (1997, 
eq. 11) and Giersz (2001, section 3.2) for clusters with masses 


around Me = 10® to IO^Mq, despite the fact that the initial 
mass-functions used and the way stellar evolution is modeled 
differ slightly. For larger clusters, we predict lifetimes which 
are smaller by a factor of 2 to 5 since eq. 10 leads to a scaling 
of the cluster lifetime with cluster mass flatter than a scaling 
with the relaxation time. 

As an application of our results for the lifetimes, we 
will have a look at how dynamical evolution depletes the 
galactic globular cluster system. Fig. 25 compares masses 
and distances of galactic globular clusters with the minimum 
mass needed such that a cluster survives for another Hubble 
time. The distances and masses are again taken from Harris 
(1996). The minimum mass was determined by first using eq. 
10 to calculate the mass of a cluster which survives for two 
Hubble-times and then eq. 12 to determine its present day 
mass. We calculated lifetimes by assuming that all clusters 
either move on circular or eccentric orbits with e = 0.5. 
For the eccentric orbits we used the fact that the ratio of 
the time-averaged galactocentric distance of a cluster in an 
e = 0.5 orbit to its apogalactic distance is equal to <Rg> 
/Ra = 0.74 if the cluster is moving through an logarithmic 
galactic potential. All present distances are divided by this 
ratio to obtain an estimate of Ra for each cluster. 

It can be seen that a large fraction of galactic globular 
clusters will not survive for another Hubble time. Depend¬ 
ing on which orbital type we assume, we find that between 
53% (circular) to 67% (eccentric orbit) of all galactic glob¬ 
ular clusters will be destroyed, most of them in the inner 
parts of the galaxy. The higher destruction rate of eccen¬ 
tric orbits is probably closer to the truth since the galactic 
globular cluster system as a whole has a radially anisotropic 
velocity distribution (Dinescu et al. 1999), implying rather 
large orbital eccentricities in the mean. 
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Figure 23. Velocity anisotropy for different mass groups as a 
function of enclosed mass for the same run as in Fig. 21. The 
differences remain small until close to final dissolution, when low- 
mass stars have more tangentially anisotropic distributions by up 
to 13 = 0.3 in the outer parts. 

Depending on which orbital type we assume, we hnd 
that currently between 7 to 10 globular clusters are de¬ 
stroyed per Gyr, which is in agreement with the value of 
5 ± 3 found by Hut & Djorgovski (1992) from observational 
estimates of the relaxation time in galactic globular clus¬ 
ters. Similarly, Fall & Zhang (2001) found by semi-analytical 
modeling that at present between 5 to 8% of all globular 
clusters are destroyed per Gyr, which corresponds to 7 to 
11 clusters. If we assume an exponential decrease with time, 
this would lead to destruction rates between 50% to 67% 
after 13 Gyr. Hence, although the assumptions about the 
underlying kinematic and galactic model are quite different 
in each paper, there is good agreement with our values. 


4.2 Evolution of the mass-function 

In this section we will apply the results of previous chapters 
to a sample of galactic globular clusters with well observed 
parameters. In recent years, the proper motions of a num¬ 
ber of clusters could be measured and tied to the reference 
frame of the Hipparcos Catalogue so that their orbits are 
relatively accurately known (see Odenkirchen et al. 1997, 
Dinescu et al. 1999 and references therein). In addition, ac¬ 
curate determinations of the mass-functions at the low-mass 
end were obtained with the HST for about a dozen globular 
clusters (Piotto & Zoccali 1999, Paresce & de Marchi 2000 
and references therein). 

Table 2 lists 10 clusters for which orbital information 
and determinations of the slope of the mass-functions are 
simultaneously available. Beside the names, we have listed 
the present-day masses, calculated from the absolute lu- 
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Figure 25. Galactocentric distance vs. mass of galactic globular 
clusters. The solid line shows the minimum mass needed such that 
a cluster survives for another Hubble time, assuming it moves in 
a circular orbit. The dashed line shows the same for clusters in 
eccentric orbits with e = 0.5 if one assumes that each cluster 
is located at the mean galactocentric distance of such an orbit 
<Rg>= 0.74:Ra- Clusters which will be destroyed in the eccentric 
case are shown by crosses, clusters which survive by dots. In the 
eccentric case, 67% of all globular clusters will be destroyed within 
the next Hubble time. 


minosities given in the compilation of Harris (1996) and 
an assumed mass-to-light ratio of 2. The perigalactic and 
apogalactic distances of the clusters in columns 3 and 4 were 
taken from Dinescu et al. (1999), except for Pal 5, where 
the data for the mass and the orbital information are from 
Odenkirchen et al. (2001, 2002). Golumn 5 lists the observed 
slope of the mass-function at the low-mass end, determined 
mostly from observations taken around the half-light radius. 

Column 6 lists estimated slopes for the global mass- 
function, determined by correcting the local mass-function 
for mass-segregation with the help of multi-mass King- 
Michie models. Since this method assumes energy parti¬ 
tion throughout the cluster, and since we have seen that 
energy equipartition is not established in the center until 
core-collapse, these corrections could be in error for some 
clusters. However, as can be seen in Table 2, the corrections 
made are usually fairly small and should therefore not sig- 
nihcantly influence our conclusions. The values in columns 
5 and 6 come from Piotto & Zoccali (1999), except for Pal 

5 and NGC 6121, where we have taken data from Grillmair 

6 Smith (2001) and Bedin et al. (2001) respectively. For 
both clusters, observations were done around the half-mass 
radius. In addition, Grillmair & Smith (2001) did not hnd 
strong evidence for mass-segregation in case of Pal 5. Hence 
we assume for both clusters that the global mass-functions 
have the same slopes as the measured ones. 

In order to make predictions for individual clusters, we 
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Table 2. Observed and predicted data for some galactic globular clusters. 


Name 

Me 

\Mq] 

Rp 

[kpc] 

Ra 

[kpc] 

Xlocal 

X global 

Foiss 

[GYR] 

Mo 

[M©] 

For 

Xth 

NGC 6121 

1.30 

■ 10® 

0.6 

5.9 

0.0 

0.0 

17.3 

7.5 

■ 10® 

0.46 

-0.32 

NGC 6254 

1.67 

■ 10® 

3.4 

4.9 

-0.1 

-0.2 

34.3 

3.8 

■ 10® 

0.33 

0.14 

NGC 6341 

3.26 

■ 10® 

1.4 

9.9 

-0.1 

0.2 

35.9 

7.3 

■ 10® 

0.32 

0.16 

NGC 6397 

7.00 

■ 10^ 

3.1 

6.3 

-0.5 

-0.35 

24.3 

2.2 

■ 10® 

0.37 

0.00 

NGC 6656 

4.29 

■ 10® 

2.9 

9.3 

-0.1 

-0.1 

62.2 

7.7 

■ 10® 

0.30 

0.25 

NGC 6752 

2.11 

■ 10® 

4.8 

5.6 

-0.3 


46.6 

4.2 

■ 10® 

0.31 

0.21 

NGC 6809 

1.79 

■ 10® 

1.9 

5.8 

-0.35 

-0.35 

28.6 

4.7 

■ 10® 

0.34 

0.08 

NGC 7078 

7.96 

■ 10® 

5.4 

10.3 

0.1 

0.4 

135.6 

1.3 

■ 10® 

0.29 

0.29 

NGC 7099 

1.59 

■ 10® 

3.0 

6.9 

0.2 

0.25 

34.7 

3.6 

■ 10® 

0.32 

0.15 

PAL 5 

5.5 

10® 

7.0 

19.0 

-0.5 

-0.5 

17.2 

3.2 

■ 10^ 

0.46 

-0.32 


will neglect the influence of a galactic disc and assume that 
clusters move through a spherical logarithmic potential. We 
then proceed similar to the previous section and calculate 
for each globular cluster the eccentricity of its orbit from 
columns 3 and 4 of Table 2 and estimate its dissolution time 
and its initial mass by using eqs. 10 and 12. These values 
can be found in columns 7 and 8 of Table 2. 

The cluster Pal 5 has the smallest lifetime in the sam¬ 
ple according to our predictions. This is also confirmed by 
observations, since Odenkirchen et al. (2001, 2002) find evi¬ 
dence that this cluster has strongly populated tidal tails and 
a very low mass, which both indicate that it must be close 
to final dissolution. 

From a comparison of the dissolution time with the clus¬ 
ter ages, the global slope of the mass-function of low-mass 
stars and the dark-matter content were finally determined 
by using eqs. 13, 15 and 16. These values can be found in the 
last two columns of Table 2. Since our simulations do not 
consider the effect of disc-shocking, our dissolution times 
are likely to be upper limits and the changes that the mass- 
functions underwent could on average be larger than esti¬ 
mated by us, especially for clusters which are close to final 
dissolution. The differences might however not be very large: 
Gnedin & Ostriker (1997) for example found that the mean 
destruction rate is increased by about 30% compared to the 
pure relaxation case if their simulations include a galactic 
disc. According to their simulations, the influence of a disc 
is strongest for a cluster near the galactic center and de¬ 
creases further outward. At Rg = 3 kpc for example, which 
corresponds roughly to the median perigalactic distance of 
the clusters in Table 2, the inclusion of a disc increases the 
average destruction rate by only 10%. Similarly, the data 
listed by Vesperini & Heggie (1997) in their Table 2 indi¬ 
cates that the influence of disc-shocking is diminishing with 
increasing distance from the galactic center and seems to 
enhance the mass-loss rate at distances of several kpc by no 
more than 10%. An inclusion of a galactic disc would there¬ 
fore probably change our results by only a small amount. 

Fig. 26 depicts the observed and global slopes of the 
stellar mass-function for the clusters in Table 2 in depen¬ 
dence of the fraction of time spent until complete cluster 
dissolution. Our simulations predict that the slope should 
decrease constantly as the clusters evolve, and there is in¬ 
deed a good correlation between the slopes and the cluster 
lifetimes in both panels. The only discrepant cluster is NGC 
6121, which is close to final dissolution but still has a flat 
mass-function with slope x = 0. Since this cluster has the or¬ 


bit with the largest eccentricity of all clusters in the sample, 
its small dissolution time is due mainly to its large orbital 
eccentricity and could be significantly larger if the true or¬ 
bit is less eccentric. We also note that Richer et al. (2002) 
have determined a mass-function slope for NGG 6121 which 
corresponds to a; = —0.25 in our notation. Using their value 
would remove most of the discrepancy with our theory. 

A Spearman rank-order correlation coefficient test 
shows that the observed and global slopes are correlated 
with the cluster lifetimes by confidence levels of 82% and 
97%. Although one should be cautious since the sample size 
is still rather small, we conclude that dynamical evolution 
has depleted the stellar mass-function in globular clusters. 

The solid line in the right panel shows our prediction for 
the global mass-function. It is for all dissolutionary stages 
within the range of observations, so it seems quite possi¬ 
ble that the galactic globular cluster system has on average 
started with a slope of a; = 0.3 at the low-mass end, similar 
to what Zoccali et al. (2000) and Kroupa (2001) have found 
for the mean mass-function of low-mass stars in the galac¬ 
tic bulge and disc. At any dissolutionary stage, the scatter 
among different clusters is considerable and seems to exceed 
the errorbars, which would hint to additional parameters 
which influence the mass-function. Next to binaries, metal- 
licity effects have long been thought of influencing the mass- 
function (McClure et al. 1986). We have therefore plotted 
the metallicity from the McMaster database next to each 
cluster in the right panel. It can be seen that all clusters 
with negative global slopes have small metallicities while 
clusters with metallicities larger than [Fe/H] = — 2 have 
slopes around x = 0.3. A second parameter which influences 
the stellar mass-function in globular clusters might therefore 
be the cluster metallicity. However, given the small number 
of data points and the fact that the mass-functions were 
deduced by different authors using various methods, this 
conclusion is still uncertain. 


5 CONCLUSIONS 

We have performed a set of large scale N-body simulations of 
multi-mass star clusters moving through an external galaxy 
with a logarithmic density profile. Our simulations included 
a realistic spectrum of stellar masses, mass-loss due to stel¬ 
lar evolution, two-body relaxation and a fully consistent 
treatment of the external tidal field. The simulations were 
carried out with the collisional N-body code NBODY4 on 
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Figure 26. Slopes of the mass-function against the fraction of time spent until complete dissolution for the clusters of Table 2. Observed 
mass-functions are shown in the left panel, global mass-functions are shown in the right panel. Errorbars are taken from Piotto & Zoccali 
(1999) and show the uncertainties of the observed slopes. The solid line in the right panel shows our theoretical prediction. There is a 
clear correlation between the local and global slope of the mass-function and the dissolutionary stage of the cluster. 


the GRAPE6 boards of Tokyo University and thanks to the 
speed of these boards we could perform simulations with 
particle numbers far greater than was possible hitherto. Our 
largest clusters have lifetimes well in excess of a Hubble time 
and become comparable to small globular clusters in the 
later stages of their evolution. 

Our results can be summarized as follows: All simulated 
clusters dissolve mainly as a result of two-body relaxation 
and external tidal truncation, with the lifetimes of clusters 
moving on the same orbit scaling proportional to the relax¬ 
ation time to a power x that is found to be in the range 0.75 
to 0.85. This is in agreement with the results of Baumgardt 
(2001) for the scaling of the lifetimes of single-mass clusters 
in tidal fields. Clusters with stronger initial concentration 
Ho show slightly larger exponents than less concentrated 
clusters. The scaling law does not change significantly if one 
goes from the circular case to clusters moving on eccentric 
orbits, which is a hint that relaxation is the dominant disso¬ 
lution mechanism for eccentric orbits as well. Utilizing our 
results for the lifetimes, we predict that about 2/3 of all 
galactic globular clusters will be destroyed within the next 
Hubble time. This value is in good agreement with recent 
literature values and argues for a strong evolution of the 
galactic globular cluster system. 

In agreement with Vesperini & Heggie (1997), Taka- 
hashi & Lee (2000) and what one would expect from the 
effect of mass-segregation, we find that low-mass stars are 
preferentially depleted from star clusters. The depletion is 
strong enough to turn an initially increasing mass-function 
into one which is decreasing towards low-mass stars. We find 
that the details of this process are nearly independent of the 
starting condition like cluster orbit, number of cluster stars 
or central concentration, and can be characterised by a sin¬ 


gle parameter, as e.g. the fraction of time remaining until 
complete cluster dissolution or the fraction of mass already 
lost from a star cluster. 

Mass segregation causes mass-functions to be steeper 
in the inner parts than near the tidal radius and might be 
difficult to correct for observationally since energy equipar- 
tition is established only very late in the cluster evolution. 
In the center it is not complete until clusters have gone into 
core-collapse while in the outer parts it is not reached until 
clusters are almost completely dissolved. This confirms an 
earlier result obtained by Giersz & Heggie (1997) and shows 
that the application of multi-mass King-Michie models to 
determine the global mass-function from observed local ones 
is not without risk. A better strategy might be to conduct 
observations at or slightly beyond the (projected) half-light 
radius, where we always find a good agreement between the 
observed slope of the mass-function and the global one, and 
make no correction for mass-segregation at all. 

The mass and number fractions of white dwarfs and 
neutron stars are increasing throughout cluster evolution 
since these stars are the most massive components and are 
therefore least likely to escape. For globular clusters near 
the end of their lifetime, like e.q. Pal 5, we predict that 
most part of the cluster mass should be in the form of com¬ 
pact remnants, unless the initial mass-function had a slope 
significantly steeper than Salpeter above m = 1 M©. White 
dwarfs and neutron stars are strongly concentrated towards 
the center and central mass-to-light ratios can reach val¬ 
ues of 10 or more at core-collapse time and in post-collapse 
evolution. This could explain the recent detection of high 
central mass-to-light ratios in post-collapse globular clusters 
without the need to assume intermediate mass black holes 
in their centers (Baumgardt et al. 2002). However, more 
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detailed simulations with a larger number of cluster stars 
and containing a significant fraction of primordial binaries, 
possible for example with Monte-Carlo codes (Giersz 2001, 
Spurzem & Giersz 2002), might be necessary to study this 
question in detail. 

Clusters remain isotropic in their centers and near the 
half-mass radius and have mildly tangential anisotropies in 
their outer parts. The amount of tangential anisotropy de¬ 
pends on the orbit of the cluster and is largest for clusters 
moving in circular orbits. High-mass stars have slightly less 
tangentially anisotropic velocity distributions in the outer 
cluster parts although the differences remain small until 
close to the end. In the outer parts, the simulated clusters 
also show a small amount of rotation since stars orbiting in 
the same sense as the cluster moves around the galaxy are 
more easily lost from the cluster than stars in retrograde 
motion. 

For a sample of galactic globular clusters with well ob¬ 
served parameters, we find a correlation between the ob¬ 
served slope of the mass-function and the lifetimes predicted 
by us in the sense that clusters close to dissolution are 
depleted in low-mass stars. There is also good agreement 
between our prediction for the slope of the mass-function 
at any evolutionary stage and the mean observed slope. 
It therefore seems conceivable that galactic globular clus¬ 
ters started on average with slopes around x = 0.3 for 
m < 0.7 Mq, which is similar to what was found by Zoc- 
cali et al. (2000) and Kroupa (2001) for low-mass stars in 
the galactic bulge and disc. This would indicate that no 
dramatic change has happened to the mass-function of low- 
mass stars within the last 10 to 12 Gyrs and that the mass- 
function does not depend much on the density of the envi¬ 
ronment where stars are formed. There are hints for addi¬ 
tional parameters which influence the mass-function, as e.g. 
the cluster metallicity. However, given the fact that the num¬ 
ber of globular clusters with well observed mass-functions is 
still rather small, this conclusion is uncertain. 
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